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Abstract 

This paper considers the class of Levy processes that can be written as a 
Brownian motion time changed by an independent Levy subordinator. Exam- 
ples in this class include the variance gamma model, the normal inverse Gaus- 
sian model, and other processes popular in financial modeling. The question 
addressed is the precise relation between the standard first passage time and 
an alternative notion, which we call first passage of the second kind, as sug- 
gested by Q and others. We are able to prove that standard first passage time 
is the almost sure Umit of iterations of first passage of the second kind. Many 
different problems arising in financial mathematics are posed as first passage 
problems, and motivated by this fact, we are lead to consider the implications 
of the approximation scheme for fast numerical methods for computing first 
passage. We find that the generic form of the iteration can be competitive 
with other numerical techniques. In the particular case of the VG model, the 
scheme can be further refined to give very fast algorithms. 

Key words: Brownian motion, first passage, time change. Levy subordinators, 
stopping times, financial models. 
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1 Introduction 



First passage problems are a classic aspect of stochastic processes that arise in many 
areas of application. In mathematical finance, for example, first passage problems 
lie at the heart of such issues as credit risk modelling, pricing barrier options, and 
the optimal exercise of american options. If Xt is any process with initial value 
Xq = xq, the first passage time to a lower level b is defined to be the stopping time 

tlixo) =mf{t>0\Xt<b} (1) 

The distributional properties of t* are well studied when the underlying process X 
is a diffusion, but when X has jumps, for example, for Levy processes, much less 
is known. 

Explicit formulas for first passage are known only for certain special Levy pro- 
cesses, specifically the Kou-Wang jump diffusion model [9| and its generalization 
to phase type processes [[3]|. In the setting of general Levy processes, one analyt- 
ical approach to first passage is the Wiener- Hopf method, described in [IJ. This 
method can for example compute first passage problems for processes with one- 
sided jumps. A second general approach to first passage is to solve the Fokker- 
Planck equation for the probability density of Xt conditioned on the set {t* > t}. 
For Levy processes this amounts to solving a certain linear partial integral differ- 
ential equation (PIDE) with Dirichlet boundary conditions. The FIDE approach to 
first passage in the general setting has apparently not been widely implemented: it 
seems that when faced with difficult first passage problems, practitioners often fall 
back on Monte Carlo methods. 

Our purpose here is to present a new approach to first passage problems appli- 
cable whenever the underlying Levy process can be realized as a Levy subordinated 
Brownian motion (LSBM), that is, whenever X can be constructed as o T where 

is a standard drifting Brownian motion and T is a non-decreasing Levy process 
independent of W. The class of Levy processes that are realizable as LSBMs is 
broad enough to include most of the Levy processes that have so far been used in 
finance, such as the Kou-Wang model, the variance gamma (VG) model, and the 
normal inverse Gaussian (NIG) model. 

The basis for our approach is that for processes that are realizable as time 
changes of Brownian motions, there is an alternative notion that is also relevant, 
namely, the first time the time change exceeds the first passage time of the Brown- 
ian motion. This notion, called first passage of the second kind in [|71, shares some 
characteristics with the usual first passage time and can be applied in a similar way. 
The usefulness of this new concept is that it can be computed efficiently in many 
cases where the usual first passage time cannot. 

In the present paper, we study first passage for LSBMs and show how the first 
passage of the second kind is the first of a sequence of stopping times that converges 
almost surely to the first passage time. Expressed differently, first passage can be 
viewed as a stochastic sum of first passage times of the second kind. This sequence 
leads to a convergent and computable expansion for the first passage probability 
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distribution function p* in terms of a similar function p\ that describes the first pas- 
sage distribution of the second kind. The outline of the paper is as follows. In 
Section 2, we define the objects needed to understand first passage time, and prove 
the expansion formula for first passage. In Section 3, we demonstrate the usefulness 
of this expansion by proving several explicit two dimensional integral formulas for 
p\, the first passage distribution of the second kind. Section 4 provides two proofs 
of the convergence of the expansion. The first proof is a proof of convergence in 
distribution, the second is in the pathwise (almost sure) sense. Section 5 focusses 
on the special case of the variance gamma (VG) model. In this important example, 
the formula for p\ is reduced to a one-dimensional integral (involving the exponen- 
tial integral function),. In Section 6, the expansion of the function p* is studied 
numerically, and found to be numerically stable and efficient. 

2 First Passage for LSBMs 

Let Xt be a general Levy process with initial value Xq = Xq and characteristics 
(6, c, u)h with respect to a truncation function h{x) (see [|3 or [|8l). This means X 
is an infinitely divisible process with identical independent increments and cadlag 
paths almost surely. 6, c > are real numbers and is a sigma finite measure on 
M \ that integrates the function 1 A a;^. By the Levy-Khintchine formula, the 
log-characteristic function of X^ is 

logE[e^"^*] = iuh - + I (e*"^ - 1 - xh{x)) v{dx). (2) 

In what follows we will find it convenient to focus on the Laplace exponent of X: 

iPx{u) :=logE[e-"^^] 

For simplicity of exposition, we specialize slightly by assuming that v is continuous 
with respect to Lebesgue measure u{dx) = u{x)dx, and integrates 1 A allowing 
us to take h{x) = 0. In this setting, the Markov generator of the process Xt applied 
to any sufficiently smooth function f{x) is 

[Cf]{x) = bdj + ^-dlJ + / (/(x + y)- fix)) u{y)dy (3) 
^ Jr\o 

Definition 1. For any 6 G M, the random variable tl = tKxo) := inf{t \ Xt < b} 
is called the first passage time for level b. When 6 = 0, we drop the subscript and 
t* := inf {t I X( < 0} is called simply the first passage time of X. 

Remarks 2. 1 . Since distributions of the increments of X are invariant under 
time and state space shifts, we can reduce computations of t^(xo) to compu- 
tations of t*{xQ — b). 
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2. A general Levy process is a mixture of a continuous Brownian motion with 
drift and a pure jump process, and t* is the minimum of a predictable stopping 
time (coming from the diffusive part) and a totally inaccessible stopping time 
(coming from the down jumps). Only if supp(z/) C M+ is t* predictable. If 
c = or if 1/ has an infinite activity of down jumps (i.e. if z/(]R_) = oo), then 
t* is totally inaccessible. 

3. When Xt* —Xt*_ ^ 0, we say that X jumps across 0, and define the overshoot 
to be Xt*. 

The central object of study in this paper is the joint distribution of t* and the 
overshoot Xt*, in particular the joint probability density function 

p*ixo; s,xi) = E,^[6it* - s)5iXt* - x^)] (4) 

The marginal density of t* is p*{xo; s) = P*ixo; s, Xi) dxi. 

In the introduction, we noted that general results on first passage for Levy pro- 
cesses, in particular results on the functions p*, are difficult to come by. For this 
reason, we now focus on the special class of Levy processes that can be expressed 
as a drifting Brownian motion subjected to a time change by an independent Levy 
subordinator. Such Levy subordinated Brownian motions (LSBM) have been stud- 
ied in a general review by Q and more specifically in Q. The general LSBM is 
constructed as follows: 

1. For an initial value xq > and drift /3, let Wt = Xq + Wt + PT be a drifting 
BM; 

2. For a Levy characteristic triple (6, 0,/i) with 6 > and supp(yu) C let 
the time change process Tt be the associated nondecreasing Levy process (a 
subordinator), taken to be independent of W; 

3. The time changed process Xt = is defined to be a LSBM. 

So constructed, it is known that Xt is itself a Levy process: fSl provide a char- 
acterization of which Levy processes are LSBMs. It was observed in that for 
any LSBM Xt, one can define an alternative notion of first passage time, which we 
denote here by i. 

Definition 3. For any LSBM Xt = we define the first passage time ofW to be 
T* = T*{xo) = inf{T ■.xo + Wt + PT< 0}. Note T*{xo) = when xq < 0. The 
first passage time of the second kind of Xt is defined to t = t(xo) = inf{t : Tt > 
T*ixo)} 

This definition of i, and its relation to t*, is illustrated in Figure [2j 
We now show that i(xo) := tl{xo) is the first of a sequence of approximations 
{t*(xo)}j=i,2,... to the stopping time t*. The construction of t*(xo,7r) is pathwise: 
We introduce the second argument tt G which denotes a sample path, that is. 
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Figure 1: Three trajectories of Xt and a sample path of time change r< 



a pair (cj, r) where u is a continuous drifting Brownian path W and r is a cadlag 
sample time change path T. Thus vr : (S, s) {^(S), t{s))s,s>o- The natural "big 
filtration" {J-'t)t>o for time-changed Brownian motion has J-'t = (j{uj{S), t{s), S < 
r(t), s < t}. For any t > there is a natural "time translation" operation on paths 
pt : (uj, t) {uj', t') where uj'{S) = uj{S + r(t)), r'(s) = r(s + t) - rif). 

The construction of {t*(xo, vr)}j=i 2,... for a given sample path vr is as follows. 
Inductively, for i > 2, we define the time of the i-th excursion overjump 

t*{xo, vr) = inf{t > tUixo, ir) : - Tt^_^ > r*(X,._^, vr')} (5) 

where vr' = pt*_^(vr) denotes a time shifted sample path. Note that t*(xo,vr) = 
t*_]^(xo, vr) if and only if Xt*__^(xo,Tr) < or t*_^(xo, vr) = 00. At any excursion over- 
jump event t*, the time interval which covers the event has left and right endpoints 
T-" = Tj._ and T^^ = Tt*. The joint distribution of t*{xo), Xf-i^xo) can be written 

p*(xo; s,x) = E^Mt: - s)SiXt: - x)] (6) 

The definition of this sequence of stopping times is summarized by the pathwise 
equation 

^i(a;o,vr) = t*_i(a;o,7r)l|x,._^<o}+(iI(^t-_,,vr') +d(xo,vr)j l{x,*_^>o}, i > 2 

(7) 

where n' = pt*_^{n). The identical increments property of the LSBM implies the 
joint probability densities satisfy the recursive relation 

00 s 

p*{xo;s,x) = pI{xo;s,x)I{x <0)+ j dy j du p*_^{xo;u,y)pl{y; s-u, x), i>2 



(8) 
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Similarly, the PDF of the first passage time t* satisfies the relation 



p*{xo;s)= J pl{xo;s,x)dx + j dy j du p*_^{xo;u)pl{y; s - u,x), i > 2. 

-oo 

(9) 

Examples of LSBMs: We note here three classes of Levy processes that can be 
written as LSBMs and have been used extensively in financial modeling. 

1. The exponential model with parameters (a, b, c) arises by taking Tt to be the 
increasing process with drift 6 > and jump measure jj,{z) = ce'""^, c, a > 
on (0, oo). The Laplace exponent of T is 

iIjt{u) := -logE[e~"'^i] = bu + uc/{a + u). 

The resulting time-changed process Xt := has triple {/3b, b, p) with 

n{y\ = ^ ^-(v//32+2a+/3)fa)+-(V/32+2a-/3)fa)- 

where (?/)+ = max(0, y), (?/)~ = [—y)^. This forms a four dimensional 
subclass of the six-dimensional family of exponential jump diffusions studied 
byia. 

2. The variance gamma (VG) model flO^ arises by taking Tt to be a gamma 
process with drift defined by the characteristic triple (6, 0,/i) with 6 > 
(usually b is taken to be 0) and jump measure [i{z) = {i^z)~^ exp(— z/i^), u > 
on (0, oo). The Laplace exponent of Tt,t = 1 is 

iPt{u) := -logE[e~"^i] = bu + ^log(l + uu). 
The resulting time-changed process has triple (pb, b, p) with 



p{y) = ^ ( ~ Y ~ + ^^1^1 



3. The normal inverse Gaussian model (NIG) with parameters /3, 7 [!45 arises 
when Tt is the first passage time for a second independent Brownian motion 
with drift > to exceed the level 7t. Then 

and the resulting time-changed process has Laplace exponent 



iJx{u) = xfi + 7[/3 + ^J P'' - + 2pu\. 
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3 Computing First Passage of the Second Kind 



We have just seen that first passage for LSBMs admits an expansion as a sum of 
first passage times of the second kind. In this section, we show that this expansion 
can be useful, by proving several equivalent integral formulas for computing the 
structure function pKxq; s, Xi) for general LSBMs. While the equivalence of these 
formulas can be demonstrated analytically, their numerical implementations will 
perform differently: which formula will be superior in practice is not a priori clear, 
but will likely depend on the range of parameters involved. For a complete picture, 
we provide independent proofs of the two given formulas. 

Theorem 1. Let ip{u) be the Laplace exponent ofTi, and let W have drift P ^ 0. 
Then 

pI{xo; s,xi) 

47r2 J J tizi-Z2) ^(3^ - 2iZ2 

2e/3(^.-^o)^^ r r ^^^^^^ cos \x, \ A:, sin x^h (^2 + ^.y,^ 



71-2 yy{M+)2 kf - ki 



TX 



sin \xi\k sinxofc e-'^^^'^ ^/'^'^^{{k'^ + (3'^)/2)dk (11) 



Here PV denotes that the principal value contour is taken. 

Remark 4. The equivalence of these two formulas can be demonstrated directly 
by performing a change of variables kj = i^(3'^ — 22-2^, j = 1,2, followed by 
a deformation of the contours. Justification of the contour deformation (from the 
branch of a left-right symmetric hyperbola in the upper half /c^ -plane to the real axis) 
depends on the decay of the integrand, and the computation of certain residues. 



3.1 First proof of Theorem [T] 

For a fixed level h > 0, the first passage time and the overshoot of the process Tt 
above the level h are defined to be i(h) = inf {t > \ Tt > h} and 5{h) = Tk^h) " ^• 
The Pecherskii-Rogozin identity [.11,1 applied to the nondecreasing process T says 
that 



-z\h 



E 



-Z2&{h) — z-ii{h) 



dh 



Zl - Z2 



[z:i + ip{zi)) V 



Inversion of the Laplace transform in the above equation then leads to 



E 



-Z25{h)~zzi{h) 



1 

2^ 



i){iZi) - %1){Z2) 

izi - Z2 



{z, + iJi^z^))-' e'^^Uz^ (12) 
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The first passage time of tlie BM witli drift is defined as T* = T*(xo) = 
inf{T > \ xq + Wt + I3T < 0}. Next, we need to find the joint Laplace 
transform of tl = inf{t | > T*} = i{T*) and the overshoot S* ^ S{T*). Since 
Tt is independent of Wt we find that 



^ ^^-Z2s*-Z3t^ ^ E E e-^25{T*)-z3i{T*) I 
1 f ip{izi) - ij{z2 



(13) 



27r 
1 

2^ 



ijjizi) - -0(2:2) 

iZi — Z2 



{z3 + ,ljiizi)y'E[e''^^'']dzi 
{z, + ^{^z,)re'<^^^^^)dz„ 



where in the last equality we have used the following well-known result for the 
characteristic function of the first passage time of BM with drift: 



Next we use the Fourier transform of the PDF of the BM with drift in time 
variable to obtain 



E[S{Wt-x,)] 



e 2t 



-iz2t 



-\xi\^0^-1izi 



'l-Kt 27r j " _ 2iz2 

R 

Thus, using the fact that W is independent of and 6* we obtain 



-dzo 



(14) 



E 



27r 



E 



E 



-1x11^/32-2122 
Z3tl-tz,5'^ d22 



47r2 Jj i{zi-Z2) y;^^^^2?^ 

R2 ^ 

Now, the statement of the Theorem follows after one additional Fourier inversion: 



pI{xo; s,xi 
1 

2^ 



E 



E 



S{tl-s) (Ws*-xi^ 
e-'''^*^5(Ws^ -xi^ dzs 



47r2 J J i{zi-Z2) Jp^ - 2iz2 

R2 

where we have also used the following Fourier integral 



-xoV /32 -2izi -\xi\^J /32 -2iz2 



2x7 i 



27r J ^2:3 + ip{izi) 



dz. = e-*^(^"il 



n 



3.2 Second proof of Theorem [I] 

The strategy of the proof is to compute 

I{u) = i?o,xo[l{s<t*<s+u}5(-^s+« — Xi)] 



(15) 



and then take the limit of I{u)/u as m — » 0+. The key idea is to note that Xg+u = 
Xs_ + W^, where W', T' are copies of W , T, independent of the filtration jFs_. We 
can then perform the above expectation via an intermediate conditioning on jFs_: 



{s<t'[<s+u} 

= l{s<q}E[6i£ + W^^-Xi)l{u>t'-t}\X,_ 



(16) 
(17) 



To evaluate the expectations that arise, we will need the second and third of the 
following results that were stated and proved in [|71 : 



Lemma 2. 1. For any s > 



Eo,x[i-{s<ti}S{Xs-y)] = l{jy>0}- 



^J3(y-x) 

2ti 



2. For any s > and e G 

E{)^x['^{s>tl}5{Xs — 



2ti 



(18) 



Jz{x+\y\)^-s^iiz^+f3^)/2)^^_ ^^^^ 



3. For any k in the upper half plane, 



First, using ( 19), we find 



E[5{Xs+u — Xi)l{s<tl<s+u}\Xs 

2n 



L{s<q}- 



k — z k + z 

(20) 



^ i] (21) 

^^gife(^+|xi|)g-«^((fc2+/32)/2)_ ^22) 



When we paste this expression into the final expectation over Xg^ we can use Fubini 
to interchange the expection and integral providing we choose e > 0. Then we find 

JR+ie 



We can now use ([20]) from Lemma [2] obtain 

/d{xi-xo) 



(27r) 



^ik\xi\+izxo 



\-ie)xl 



k — z k + z 



(24) 



Noting that /(O) = and taking lim„_^o I{u)/u now gives 



P*i{xq] S,Xi) 

g/3(xi-a;o) 

27r2 



(25) 



(IR+ie)xI 



K"' — Z^ 



Here the arbitrary parameter e > can be seen to ensure the correct prescription for 
dealing with the pole at k"^ = z"^. 

Finally, the complex integration in ([26]) can be expressed in the following man- 
ifestly real form: 



2qP(xi-xo) 

g/3(xi-a;o) 
TT 



-PV 



dkdz 



. + \2 



z COS sin xqz 
k"^ — z"^ 



sm 



\xi\z sinxoz e-^'^«"'+^')/2)^^^2 ^ (^^)i2)dz. (26) 



involving a principle value integral plus explicit half residue terms for the poles 

k = ±z. □ 



4 The iteration scheme and its convergence 



The next theorem shows that ([8]) can be used to compute p*{xq] s,x). We define a 
suitable norm for functions /(xq; u, x): 



oo = sup 

xo>0 




|/(xo; u, x)\dudx 



.0 



(27) 



Theorem 3. The sequence (p* )n>i converges exponentially in the L°° norm. 
Proof. First we find from ([8]) 



oo s 



Pn+lK^O] S,Xi) -Pn[Xo; S,Xi^ 

thus 



dy / du pI{xo; s-u,y) [pKy; u, xi) - p^^iiy; u, xi)] , 







Pn+1 -Pnlloo < C\\pI-pI_^\ 
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where 



C = sup 

a;o>0 



oo oo 

J J Pi{xo;u,x)dudx 



,0 



The proof is based on the probabilistic interpretation of the constant C: by 
definition pl{xo; u, x) is the joint density of tl and Xq, thus we obtain: 

C = sup P {tl < +00, Xq>0\Xo^ Xq) . 
xo>0 

Next, using the fact that Wt* = (T* is the first passage time of Xt and X is a 
continuous process) and the strong Markov property of the Brownian motion we 
find that 

C = P(tl< +00, Wt^, - Wt* > I Wo = Xq^ 
= P(tl< +00, Ws* + > I Wo = 0) , 

where the Brownian motion Wt is independent of Tj and 6* = 6*{xo) = Tq — T* is 
the overshoot of the time change above T*. 
Thus we need to prove that 

C = sup P {tl < +00, W5. + P6* > I Wo = 0) < 1, (28) 

a;o>0 

where tl — tl{xo), 5* — S*{xo) and the Brownian motion W is independent of tl 
and S*. 

First we will consider the case when < 0. In this case we obtain 

P {tl < +00, Wa* + > I Wo = 0) < P {Ws* + > | Wo = 0) 

00 00 

= J P{Wt + (3t > \Wo ^0)P{5* e dt) < J^P{5* edt)^^, 


where we have used the fact that Wt is independent of the overshoot 5* and that 
P{Wt + /5t > I Wq = 0) < I for any t and any /? < 0. Thus in the case when the 
drift P is negative we obtain an estimate C < |. 

The case when the drift P is positive is more complicated. We can not use the 
same techniques as before, since the bound P{Wt + pt > \ Wq = 0) < | is no 
longer true: in fact P{Wt + pt > \ Wq = 0) monotonically increases to 1 as 
t ^ 00. 

First we will consider the case when xq is bounded away from 0: a^o > c > 0. 
Then xo + Wt + pt has a positive probability of escaping to +00 and never crossing 
the barrier at 0, thus 

P {tl < +00, Ws* +p6* >0\Wo = 0) < P (tlixo) < +00) 

< P(ri(c) < +00) = l-ei(c). 
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Now we need to consider the case when xq — >• 0+. The proof in this case is 
based on the following sequence of inequalities: 

P {tl < +00, Ws* +(36* >0\Wo = 0) (29) 
< P {Ws* +(36* >0\Wo = 0) 
= 1-P {Ws-^ + (36* <0\Wo = 0) 



1 - 



j P{Wt + (3t<Q\Wo = 0)P{6* G dt) 


r 

j P{Wt + pt<0\Wo = 0)P((5* G dt) 



< 1 

< 1- J P{Wr + pT<0\Wo = 0)P{6* G dt) 



= 1- P{Wr + (3r <0\Wo = 0)P{6* <t), 

where r is any positive number and the last inequality is true since P{Wt + (3t < 
I Wq = 0) is a decreasing function of t. 

Since Xq 0+, we also have T*{xq) 0+ with probability 1. Since 6* is 
the overshoot of T*, and T* 0"*" as Xq 0+, we see that the distribution of the 
overshoot 6*{xq) converges either to the distribution of the jumps of Tt if the time 
change process Tt is a compound Poisson process or to the Dirac delta distribution 
at if Tt has infinite activity of jumps. Therefore in the case when Tt is a compound 
Poisson process with the jump measure u{dx) we choose r such that z/([0, r]) > 0, 
and if Tt has infinite activity of jumps we can take any r > 0. Then we obtain 
lim P{6*{xq) < t) = ^, where ^ = z/([0,r]) in the case of compound Poisson 

process, and ,^ = 1 in the case of the process with infinite activity of jumps. Using 
(|29l) we find that as xq 0+ 



P (tlixo) < +00, W5*+(36* >0\Wo = 0) < 1 - P{Wr + (3t <0\Wo = 0)^ 

< 1-62. 

To summarize, we have proved that the function 

P(xo) = P (tlixo) < +00, Ws* +(36* >0\Wo = 0), 
satisfies the following properties: 

• for any c > there exists ei = ei(c) > such that P{xo) < 1 — ei(c) for all 

Xq > C 

• there exists €2 > such that lim P{xq) < 1 — 62 ■ 
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Therefore we conclude that there exists e > 0, such that P{xo) < 1 — e for all 
xq > 0, thus C < 1 — e. This ends the proof in the second case /3 > 0. □ 



For a complementary point of view, the next result shows that the sequence 
{t*{xo)i>i converges pathwise. 

Theorem 4. For any TCBM with Levy subordinator Tt and Brownian motion with 
drift P the sequence of stopping times (t*(a;o)i>i converges a.s to t*. 

Proof: If t* = oo, then certainly t* oo, so we suppose t* < oo. In this case, 
if t* = t*_^_'^ for some i the sequence converges, and thus the only interesting case 
to analyze is if t* 7^ t* for all i < 00. Then we have < < ■ ■ ■ < t* < 
. . . . Correspondingly, we have an infinite sequence of excursion overjump intervals 
which do not overlap: let their endpoints be T*_ := Tt*- < T*j_ := Tq. The 
following observations lead to the conclusion: 

1 . by monotonicity and boundness of the sequences {T*_ ) and ) , limj^oo T*_ = 
lim^^oo T*^ = Too exists; 

2. Xq + Wt^ + PToo = by the continuity of Brownian motion; 

3. limj_,oo t* = too exists, and too 

4. Jump times are totally inaccessible, so there is no time jump at time too almost 
surely, hence Tt^ = Too', 

5. Therefore X^^ = and so too >t*: hence too = t*. 

□ 



5 The Variance Gamma model 

The Variance Gamma (VG) process described in Section 2 is the LSBM where the 
time change process Tt is the Levy process with jump measure fi(z) = (i^z)^^ exp{—z/i') 
on (0, 00) and the Laplace exponent iPt{u) = I log(l + uu). In this section we take 
the parameter 6 = 0. This model has been widely used for option pricing where it 
has been found to provide a better fit to market data than the Black-Scholes model, 
while preserving a degree of analytical tractability. The main result in this section 



reduces the 2D integral representation (10) for pl{xo; s,xi) to a ID integral and 



leads to greatly simplified numerical computations. 



The authors would like to thank Professor Martin Barlow for providing ths proof. 
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Theorem 5. Define a = W ^ + /S^. Then 



-2iz 



X 



e-\^^\V^^'Ei (^-\xi\(^a- ^p^-2iz^^ 
where Ei{x) is the exponential integral function (see f6\]). 



(30) 



Proof. Consider the function I{zi) which represents the outer integral in ( 10 1 

1 f log(l + iuzi) - log(l + iuz2) e-l^ilV^^^2^ 



dzo 



271 J i{zi - Z2) J3'i - 2iz2 

R 

First we perform the change of variables u = %\/ — 2iz2 and obtain 



1 /• log(l + |(^^ + /3^))-log(l + ^^^i) ^.|,,|. 



TT 



m2 + /?2 _ 2iz^ 



(31) 



where the contour L obtained from M under map Z2 u = IS trans- 

formed into contour M (this is justified since the integrand is an analytic function in 
this region for any zi). To finish the proof we separate the logarithms 



log ^1 + ^(m^ + j - log(l + iuzi) = \og{u + ia) + log(n - ia) - log ^ 



2izi 



use the partial fractions decomposition 



1 



+ (3^-2izi 2i^f3^ - 2tzi 



u 



i^/(3^ - 2izi u + i^/f3^ -2izi 



and we obtain the six integrals, which can be computed by shifting the contours 
of integration and using the following Fourier transform formulas (see Gradshteyn- 
Ryjik...) 



log(l + |je-^ 



log 



1 e 



txy 



dy 

y 

dy 



y 



-2mEi{-bx), b>0 
-2mEi{bx), 6 > 



i{b+e 



□ 
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Remark: Using the change of variables u = i\/ — 2iz and simplifying the 
expression we can obtain a simpler formula for Pi(xo; s, xi): 



Pl{xo;s,Xi) = e^(^'^-^o)^^^— / (a^ + y^) ^ sin (xoy) e'^'^'^^Ei {- \ xi \ {a + ly)) dy. 



2m 



Applying the Plancherel formula to the above expression gives us the following 
representation for p\ : 



pl{xo;s,x,) = e/3(-i-o)-a(.o+|xi|)^^^ ^^tI^ (32) 

^/^Tly{xo+ I xi I) r + 1) 



+ ./^e^(^i-^o)-°I^^I^^^^V^ I ut-'^K^^^{au)f{xo,x^■u)du, 



s 1 

r(e) ./ 





where 



f{xo,xi;u) = — — r - sign (u - xo) I r— r-2 

U + Xo+ \ Xi \ \u — Xo\+ \ Xi \ Xo+ I Xi I 

The above expression is useful for computations when s is small. In particular, 
when s = we find 

^(3{xi-xo)-a{xo + \xi\) 

Pi(xo;0,xi) = — — TT — . (33) 

u{xo+ I xi I) 

6 Numerical implementation for VG model 

The algorithm for computing the functions p*{xo; s, x) and p*{xo; s) can be sum- 
marized as follows: 

1 . Choose the discretization step sizes 5^, St and discretization intervals [—X, X] , 
[0,T]. The grid points are U = i6t, 1 < i < Ns and xj = {j + l/2)5x, 
-N, < J < iV,. 



2. Compute the 3D array p\{xi]tj,Xk)- For j > use equation (30) and for 



j = use explicit formula (33 1. 



3. Iterate equation ([8]) or (|9]). This step can be considerably accelerated if the 
convolution in w-variable is done using Fast Fourier Transform methods. We 
used the midpoint rule for integration in the y and u variables. 

Theorem |3] implies that Step 3 in the above algorithm has to be repeated only 
a few times: In practice we found that 3-4 iterations is usually enough. An im- 
portant empirical fact is that the above algorithm works quite well with just a few 
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Figure 2: The density of the first passage time for the two set of parameters. The 
circles show the "exact" result; the three solid lines show the first three approxima- 
tions. 



discretization points in the x-variable. We found that if one uses a non-linear grid 
(which places more points Xj near x = 0) then the above algorithm produces rea- 
sonable results with values of A^^. as small as 10 or 20. 

We compared our algorithm for the PDF p*{xo; s) to a finite-difference method 
that was implemented as follows. First we approximated the first passage time by 
its discrete counterpart: 

i* = i*{x) = mm{U : Xt^ < 0\Xo = x} 

where tj = iSt, < i < rit is the discretization of the interval [0, T]. The probabili- 
ties fi{x) = P{i* > ti\Xo = x) satisfy the iteration: 

fi+i{x) = l^>o / p{St,x- y)fi{y)dy,i > 1 (34) 



with /o(a;) = lx>o and can be computed numerically with the following steps: 

1. discretize the space variables x = i6x, y = jS^, < i,j < n^; 

2. compute the array of transitional probabiUties pi = p{5t, Xi), and normalize 
Po so that EiPi = 1; 



3. use the convolution (based on FFT) to iterate equation (34) times; 
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.^I • , • • : -4I • • ^ • ' ' 

2 4 6 8 10 12 2 4 6 3 10 12 



Figure 3: Vertical axis is the log^Q (||j9* — P*\\li), horizontal axis is the number of 
iterations. = 10 (left) and A^^ = 20 (right) and A^^t G {10, 25, 50, 100, 200} 



4. compute the approximation of the first passage time density p*{x,ti + 6t/2) = 
(/,+i(x)-/.(x))M. 

The big advantage of this method is that it is explicit and unconditionally stable: 
we can choose the number of discretization points in x-space and t-space indepen- 
dently. This is not true in general explicit finite difference methods, where one 
would solve the Fokker-Plank equation by discretizing the Markov generator and 
derivative in time, since 5t and have to lie in a certain subset in order for the 
methods to be stable. 

Figure |6] summarizes the numerical results for the PDF p*{xo; s) over the time 
interval [0, 5] for the VG model with the following two sets of parameters: 

Set I: xo = 0.5, p = 0.2, = 1 
Set II: Xo = 0.5, p = -0.2, u = 2. 

The number of grid points used was Nt = 50 and A^^,. = 10. The red circles corre- 
spond to the solution obtained by a high resolution finite difference FIDE method 
as described above (with nt = 1000 and = 10000), and the black lines show 
successive iterations p*{xo,t) converging to p*(xo,t). As we see, 3 iterations of 
equation ([9]) provide a visually acceptable accuracy in a running time of less than 
0.1 second (on a 2.5Ghz laptop). 

Figure [6] illustrates the convergence of our method and Table 1 shows the com- 
putation times (on the same 2.5Ghz laptop). We used Set II of parameters for the 
VG process, and the PIDE method with rit = 1000 and = 10000 to compute the 
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"exact" solution p*{xo, t). Figure [6] shows the \ogiQ of the error 

T 

\\p*-P*\\l, = J \p*ixo,t)-p*ixo,t)\dt 



on the vertical axis and the number of iterations on the horizontal axis; different 
curves correspond to different number of discretization points in t-space. The num- 
ber of discretization points in x-space is fixed at = 10 for the left picture and 
= 20 for the right picture. We see that initially the error decreases exponen- 
tially and then flattens out. The flattening indicates that our method converges to 
the wrong target (which is to be expected since there is always a discretization error 
coming from and Nt being finite). However, increasing Nt and A^^; brings us 
closer to the "target". In the table 1 we show precomputing time needed to compute 
the 3D array p*i{xi; tj,Xk) and the time needed to perform each iteration 







Nt = 10 


Nt = 25 


Nt = 50 


Nt = 100 


Nt = 200 


precomputing time 


N^ = 10 


0.0313 


0.0259 


0.0324 


0.0461 


0.0687 


each iteration 


N^ = 10 


0.0006 


0.0008 


0.0011 


0.0021 


0.0046 


precomputing time 


N^ = 20 


0.0645 


0.0612 


0.0745 


0.0868 


0.1298 


each iteration 


N^ = 20 


0.0037 


0.0045 


0.0066 


0.0120 


0.0269 



Table 1: Computation time (sec) for the new approach 



To put these results into perspective, on Figure [6] and Table 2 we present similar 
results for finite difference method. On Figure |6] we show the same logarithm of 
the error on the vertical axis, and the number of discretization points Ux on the 
horizontal axis. Different curves correspond to rit G {50, 100, 200}. The running 
time presented in Table 2 includes only the time needed to perform rit convolutions 
([34]) using the FFT. As we see, the finite difference method is substantially slower 
than our method. 





rix = 1150 


rix = 2300 


rix = 3450 


rix = 4600 


rix = 5750 


rit = 50 


0.0756 


0.2935 


0.9582 


1.8757 


2.9967 


rit = 100 


0.1456 


0.5821 


1.9397 


3.7409 


6.0026 


rit = 200 


0.2870 


1.1478 


3.8833 


7.4768 


11.9935 



Table 2: Computation time (sec) for the finite difference approach. 



7 Conclusions 

First passage times are an important modeling tool in finance and other areas of 
applied mathematics. The main result of this paper is the theoretical connection 
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Figure 4: Vertical axis is the log^o {\\p* — P*\ horizontal axis is rix. The rit is in 
the set {50, 100, 200} 



between two distinct notions of first passage time that arise for Levy subordinated 
Brownian motions. This relation leads to a new way to compute true first passage for 
these processes that is apparently less expensive than finite difference methods for 
a given level of accuracy. Our paper opens up many avenues for further theoretical 
and numerical work. For example, the methods we describe are certainly applicable 
for a much broader class of time changed Brownian motions and time changed 
diffusions. Finally, it will be worthwhile to explore the use of the first passage of 
the second kind is a modeling alternative to the usual first passage time. 
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